rm(list = ls())

## SET YOUR FILE PATH
setwd("~/Dropbox/Immigration/Analysis/Policy Diffusion/Replication Files")

library(dplyr)
library(survminer)
library(tidyr)
library(standardize)
library(survival)
library(coxme)
library(frailtySurv)
library(frailtyEM)
library(frailtypack)
library(Hmisc)
library(haven)
library(foreign)
library(coefplot)
library(stargazer)
library(simPH)
library(coefplot)
library(xtable)
library(ggplot2)
library(texreg)
library(sp)
library(sf)
require(rgdal)
library(tmap)
library(tmaptools)
library(sp)
library(leaflet)
library(haven)
library(ggplot2)
library(shinyjs)
library(RColorBrewer)
library(ivtools)
library(jtools)

data <- read_dta("dwrap_index.dta")

data.africa <- subset(data, africa == 1)
data.1980 <- subset(data, year >= 1980)
data.1990 <- subset(data, year >= 1990)
data.indp5 <- subset(data, indp5 == 1)
data.indp10 <- subset(data, indp10 == 1)
data.law <- subset(data, numdomlaws >= 1)

set.seed(8675309)

## Set Convergence Criteria
coxph.control(eps = sqrt(.Machine$double.eps), toler.chol = .Machine$double.eps,
              iter.max = 1000, toler.inf = sqrt(.Machine$double.eps), outer.max = 1000)

coxme.control(eps = sqrt(.Machine$double.eps), toler.chol = .Machine$double.eps,
              iter.max = 1000)

## TABLE 3

res1 <- matrix(NA,35,9)

m1a <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               exc1inc2_reg_lag1 +
               std_polyarchy_lag1+
               majorannintra_contig_lag1+
               ihs_aidgdp_5ma_wins+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m1a)
res1[1:2,2] <- c(m1a$coef[1],sqrt(m1a$var[1,1]))
res1[4:5,2] <- c(m1a$coef[2],sqrt(m1a$var[2,2]))
res1[33,2] <- m1a$n
res1[34,2] <- m1a$history[[1]]$c.loglik
res1[35,2] <- AIC(m1a)

m2a <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               exc1inc2_reg_lag1 +
               std_polyarchy_lag1+
               majorannintra_contig_lag1+
               ihs_aidgdp_5ma_wins+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m2a)
res1[1:2,3] <- c(m2a$coef[1],sqrt(m2a$var[1,1]))
res1[4:5,3] <- c(m2a$coef[2],sqrt(m2a$var[2,2]))
res1[33,3] <- m2a$n
res1[34,3] <- m2a$history[[1]]$c.loglik
res1[35,3] <- AIC(m2a)

m3a <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               exc1inc2_reg_lag1 +
               std_polyarchy_lag1+
               majorannintra_contig_lag1+
               ihs_aidgdp_5ma_wins+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               reglib1last3+
               leglib1last3+
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m3a)
res1[1:2,4] <- c(m3a$coef[1],sqrt(m3a$var[1,1]))
res1[4:5,4] <- c(m3a$coef[2],sqrt(m3a$var[2,2]))
res1[10:11,4] <- c(m3a$coef[11],sqrt(m3a$var[11,11]))
res1[13:14,4] <- c(m3a$coef[12],sqrt(m3a$var[12,12]))
res1[33,4] <- m3a$n
res1[34,4] <- m3a$history[[1]]$c.loglik
res1[35,4] <- AIC(m3a)

keeps <- c("exc1inc2_reg_lag1", "ihs_aidgdp_5ma_wins", "std_ihs_fractal12_alesina", "landcontig", "ihs_gdppc_lag1", "ihs_pop_lag1", "std_polyarchy_lag1", "cwbroad_lag1", "ihs_iterate_lag1", "ihs_tradeshare_lag1", "majorannintra_contig_lag1", "numlib1", "fractalcontig", "govtfracprob", "lib1_start", "lib1_end", "liberalization1", "neggdpshock_lag1", "ccode2", "year")
df = data[keeps]

fitT.LX_kin <- coxph(formula=Surv(lib1_start, lib1_end, liberalization1)~
                       exc1inc2_reg_lag1 +
                       std_polyarchy_lag1+
                       majorannintra_contig_lag1+
                       ihs_aidgdp_5ma_wins+
                       ihs_gdppc_lag1 +
                       neggdpshock_lag1+
                       ihs_pop_lag1 +
                       cwbroad_lag1 +
                       ihs_iterate_lag1 +
                       ihs_tradeshare_lag1 +
                       
                       strata(numlib1)+
                       cluster(ccode2),
                     data = df,
                     method = c("efron"),
                     model = TRUE)
summary(fitT.LX_kin)

fitX.LZ_kin <- glm(formula=exc1inc2_reg_lag1~
                     std_ihs_fractal12_alesina+landcontig+
                     std_polyarchy_lag1+
                     majorannintra_contig_lag1+
                     ihs_aidgdp_5ma_wins+
                     ihs_gdppc_lag1 +
                     neggdpshock_lag1+
                     ihs_pop_lag1 +
                     cwbroad_lag1 +
                     ihs_iterate_lag1 +
                     ihs_tradeshare_lag1, data=df)
summary(fitX.LZ_kin)

m4a <- ivcoxph(estmethod="ts", fitX.LZ=fitX.LZ_kin, fitT.LX=fitT.LX_kin, clusterid="ccode2", data=df, ctrl = F)
summary(m4a)
res1[1:2,5] <- c(2.21664, 0.59403)
res1[4:5,5] <- c(-0.46524, 0.08300)
res1[33,5] <- c(2624)

m5a <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               exc1inc2_reg_lag1*std_polyarchy_lag1+
               majorannintra_contig_lag1+
               ihs_aidgdp_5ma_wins+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m5a)
res1[1:2,6] <- c(m5a$coef[1],sqrt(m5a$var[1,1]))
res1[4:5,6] <- c(m5a$coef[2],sqrt(m5a$var[2,2]))
res1[7:8,6] <- c(m5a$coef[103],sqrt(m5a$var[103,103]))
res1[33,6] <- m5a$n
res1[34,6] <- m5a$history[[1]]$c.loglik
res1[35,6] <- AIC(m5a)

m6a <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               exc1inc2_reg_lag1*std_polyarchy_lag1+
               majorannintra_contig_lag1+
               ihs_aidgdp_5ma_wins+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m6a)
res1[1:2,7] <- c(m6a$coef[1],sqrt(m6a$var[1,1]))
res1[4:5,7] <- c(m6a$coef[2],sqrt(m6a$var[2,2]))
res1[7:8,7] <- c(m6a$coef[170],sqrt(m6a$var[170,170]))
res1[33,7] <- m6a$n
res1[34,7] <- m6a$history[[1]]$c.loglik
res1[35,7] <- AIC(m6a)

m7a <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               exc1inc2_reg_lag1*std_polyarchy_lag1+
               majorannintra_contig_lag1+
               ihs_aidgdp_5ma_wins+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               reglib1last3+
               leglib1last3+
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m7a)
res1[1:2,8] <- c(m7a$coef[1],sqrt(m7a$var[1,1]))
res1[4:5,8] <- c(m7a$coef[2],sqrt(m7a$var[2,2]))
res1[7:8,8] <- c(m7a$coef[105],sqrt(m7a$var[105,105]))
res1[10:11,8] <- c(m7a$coef[11],sqrt(m7a$var[11,11]))
res1[13:14,8] <- c(m7a$coef[12],sqrt(m7a$var[12,12]))
res1[33,8] <- m7a$n
res1[34,8] <- m7a$history[[1]]$c.loglik
res1[35,8] <- AIC(m7a)

fitT.LX_kinint <- coxph(formula=Surv(lib1_start, lib1_end, liberalization1)~
                       exc1inc2_reg_lag1*std_polyarchy_lag1+
                       majorannintra_contig_lag1+
                       ihs_aidgdp_5ma_wins+
                       ihs_gdppc_lag1 +
                       neggdpshock_lag1+
                       ihs_pop_lag1 +
                       cwbroad_lag1 +
                       ihs_iterate_lag1 +
                       ihs_tradeshare_lag1 +
                       
                       strata(numlib1)+
                       cluster(ccode2),
                     data = df,
                     method = c("efron"),
                     model = TRUE)
summary(fitT.LX_kinint)

fitX.LZ_kinint <- glm(formula=exc1inc2_reg_lag1~
                     std_ihs_fractal12_alesina+landcontig+
                     std_polyarchy_lag1+
                     majorannintra_contig_lag1+
                     ihs_aidgdp_5ma_wins+
                     ihs_gdppc_lag1 +
                     neggdpshock_lag1+
                     ihs_pop_lag1 +
                     cwbroad_lag1 +
                     ihs_iterate_lag1 +
                     ihs_tradeshare_lag1, data=df)
summary(fitX.LZ_kinint)

m8a <- ivcoxph(estmethod="ts", fitX.LZ=fitX.LZ_kinint, fitT.LX=fitT.LX_kinint, clusterid="ccode2", data=df, ctrl = F)
summary(m8a)
res1[1:2,9] <- c(2.23004, 0.60600)
res1[4:5,9] <- c(-0.64800, 0.16951)
res1[7:8,9] <- c(0.33599, 0.21713)
res1[33,9] <- c(2624)

res1 <- round(res1,3)
res1[1,1] <- c("Elite Kin Discrimination")
res1[4,1] <- c("Democracy")
res1[7,1] <- c("Elite Kin Discrimination x Democracy")
res1[10,1] <- c("Regional Liberalization")
res1[13,1] <- c("Legal Origins Liberalization")
res1[33,1] <- c("Observations")
res1[34,1] <- c("Log-Likelihood")
res1[35,1] <- c("AIC")

latex(res1, file = "")
summary(m1a)
summary(m2a)
summary(m3a)
summary(m4a)
summary(m5a)
summary(m6a)
summary(m7a)
summary(m8a)

## FIGURE 5

SimInt1 <- coxsimInteract(m6a, b1 = "exc1inc2_reg_lag1", b2 = "std_polyarchy_lag1",
                          qi = "Marginal Effect",
                          X2 = seq(-.7, .6, by=.1),expMarg = F, nsim = 1000, ci = 0.95, spin=F, extremesDrop=T)

demintkin<-simGG(SimInt1, type = "points", alpha = .4, xlab = "\n Democracy (V-DEM)",
                 ylab = "Average Marginal Effect (AME) \n of Elite Kin Discrimination \n", rug=T, lcolour = "#000000", pcolour = "#000000", method="loess")

demintkin<- demintkin + scale_x_continuous(breaks = seq(-.7, .7, by = .2)) + scale_y_continuous(breaks = seq(0, 2.5, by = .25))

demintkin<-demintkin + geom_hline(yintercept=0, linetype="dashed", color = "black") + 
  ggtitle("Probability of Policy Liberalization")+
  theme(plot.title = element_text(hjust = 0.5))

demintkin

## FIGURE A.27

pha <- coxme(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               exc1inc2_reg_lag1 +
               std_polyarchy_lag1+
               majorannintra_contig_lag1+
               ihs_aidgdp_5ma_wins+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               
               strata(numlib1) +
               (1|ccode2),
             data = data,
             ties = c("efron"))

test.ph <- cox.zph(pha)
test.ph

plot(test.ph, var=1, xlab="Time \n",  ylab="", main="β(t) for Elite Kin Discrimination", sub="Individual Test χ2=0.106", col=c("#236ffc", "#ff3654"))
dev.copy(png,'sr_elite.png', width = 4, height = 4, units = 'in', res = 300)
dev.off()

plot(test.ph, var=2, xlab="Time \n",  ylab="", main="β(t) for Democracy", sub="Individual Test χ2=0.020", col=c("#236ffc", "#ff3654"))
dev.copy(png,'sr_dem.png', width = 4, height = 4, units = 'in', res = 300)
dev.off()

plot(test.ph, var=3, xlab="Time \n",  ylab="", main="β(t) for Civil Conflict in Neighbor", sub="Individual Test χ2=3.158", col=c("#236ffc", "#ff3654"))
dev.copy(png,'sr_cwn.png', width = 4, height = 4, units = 'in', res = 300)
dev.off()

plot(test.ph, var=4, xlab="Time \n",  ylab="", main="β(t) for DAC Aid/GDP", sub="Individual Test χ2=2.451", col=c("#236ffc", "#ff3654"))
dev.copy(png,'sr_dacaid.png', width = 4, height = 4, units = 'in', res = 300)
dev.off()

plot(test.ph, var=5, xlab="Time \n",  ylab="", main="β(t) for GDP/Capita", sub="Individual Test χ2=3.466", col=c("#236ffc", "#ff3654"))
dev.copy(png,'sr_gdppc.png', width = 4, height = 4, units = 'in', res = 300)
dev.off()

plot(test.ph, var=6, xlab="Time \n",  ylab="", main="β(t) for Negative GDP Shock", sub="Individual Test χ2=0.169", col=c("#236ffc", "#ff3654"))
dev.copy(png,'sr_shock.png', width = 4, height = 4, units = 'in', res = 300)
dev.off()

plot(test.ph, var=7, xlab="Time \n",  ylab="", main="β(t) for Population", sub="Individual Test χ2=0.284", col=c("#236ffc", "#ff3654"))
dev.copy(png,'sr_pop.png', width = 4, height = 4, units = 'in', res = 300)
dev.off()

plot(test.ph, var=8, xlab="Time \n",  ylab="", main="β(t) for Civil Conflict in Policymaker", sub="Individual Test χ2=0.093", col=c("#236ffc", "#ff3654"))
dev.copy(png,'sr_cw.png', width = 4, height = 4, units = 'in', res = 300)
dev.off()

plot(test.ph, var=9, xlab="Time \n",  ylab="", main="β(t) for Transnational Terrorism", sub="Individual Test χ2=0.121", col=c("#236ffc", "#ff3654"))
dev.copy(png,'sr_terror.png', width = 4, height = 4, units = 'in', res = 300)
dev.off()

plot(test.ph, var=10, xlab="Time \n",  ylab="", main="β(t) for Trade-to-GDP Ratio", sub="Individual Test χ2=3.478", col=c("#236ffc", "#ff3654"))
dev.copy(png,'sr_trade.png', width = 4, height = 4, units = 'in', res = 300)
dev.off()

## FIGURE 6

data$kin <- NULL
data$kin <- data$exc1inc2_reg_lag1

lib50 <- coxph(formula = Surv(lib5_start, lib5_end, liberalization_50) ~
                 kin +
                 std_polyarchy_lag1+
                 majorannintra_contig_lag1+
                 ihs_aidgdp_5ma_wins+
                 ihs_gdppc_lag1 +
                 neggdpshock_lag1+
                 ihs_pop_lag1 +
                 cwbroad_lag1 +
                 ihs_iterate_lag1 +
                 ihs_tradeshare_lag1 +
                 
                 cluster(ccode2) +
                 strata(numlib5) +
                 frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                 frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
               data = data,
               method = c("efron"),
               model = TRUE)

summary(lib50)

lib60 <- coxph(formula = Surv(lib6_start, lib6_end, liberalization_60) ~
                 kin +
                 std_polyarchy_lag1+
                 majorannintra_contig_lag1+
                 ihs_aidgdp_5ma_wins+
                 ihs_gdppc_lag1 +
                 neggdpshock_lag1+
                 ihs_pop_lag1 +
                 cwbroad_lag1 +
                 ihs_iterate_lag1 +
                 ihs_tradeshare_lag1 +
                 
                 cluster(ccode2) +
                 strata(numlib6) +
                 frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                 frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
               data = data,
               method = c("efron"),
               model = TRUE)

summary(lib60)

lib70 <- coxph(formula = Surv(lib7_start, lib7_end, liberalization_70) ~
                 kin +
                 std_polyarchy_lag1+
                 majorannintra_contig_lag1+
                 ihs_aidgdp_5ma_wins+
                 ihs_gdppc_lag1 +
                 neggdpshock_lag1+
                 ihs_pop_lag1 +
                 cwbroad_lag1 +
                 ihs_iterate_lag1 +
                 ihs_tradeshare_lag1 +
                 
                 cluster(ccode2) +
                 strata(numlib7) +
                 frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                 frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
               data = data,
               method = c("efron"),
               model = TRUE)

summary(lib70)

lib80 <- coxph(formula = Surv(lib80_start, lib80_end, liberalization_80) ~
                 kin +
                 std_polyarchy_lag1+
                 majorannintra_contig_lag1+
                 ihs_aidgdp_5ma_wins+
                 ihs_gdppc_lag1 +
                 neggdpshock_lag1+
                 ihs_pop_lag1 +
                 cwbroad_lag1 +
                 ihs_iterate_lag1 +
                 ihs_tradeshare_lag1 +
                 
                 cluster(ccode2) +
                 strata(numlib80) +
                 frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                 frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
               data = data,
               method = c("efron"),
               model = TRUE)

summary(lib80)

lib90 <- coxph(formula = Surv(lib90_start, lib90_end, liberalization_90) ~
                 kin +
                 std_polyarchy_lag1+
                 majorannintra_contig_lag1+
                 ihs_aidgdp_5ma_wins+
                 ihs_gdppc_lag1 +
                 neggdpshock_lag1+
                 ihs_pop_lag1 +
                 cwbroad_lag1 +
                 ihs_iterate_lag1 +
                 ihs_tradeshare_lag1 +
                 
                 cluster(ccode2) +
                 strata(numlib90) +
                 frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                 frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
               data = data,
               method = c("efron"),
               model = TRUE)

summary(lib90)

ewlib1 <- coxph(formula = Surv(ewlib1_start, ewlib1_end, ew_liberalization1) ~
                  kin +
                  std_polyarchy_lag1+
                  majorannintra_contig_lag1+
                  ihs_aidgdp_5ma_wins+
                  ihs_gdppc_lag1 +
                  neggdpshock_lag1+
                  ihs_pop_lag1 +
                  cwbroad_lag1 +
                  ihs_iterate_lag1 +
                  ihs_tradeshare_lag1 +
                  
                  cluster(ccode2) +
                  strata(ewnumlib1) +
                  frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                  frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
                data = data,
                method = c("efron"),
                model = TRUE)

summary(ewlib1)

lib1 <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
                kin +
                std_polyarchy_lag1+
                majorannintra_contig_lag1+
                ihs_aidgdp_5ma_wins+
                ihs_gdppc_lag1 +
                neggdpshock_lag1+
                ihs_pop_lag1 +
                cwbroad_lag1 +
                ihs_iterate_lag1 +
                ihs_tradeshare_lag1 +
                
                cluster(ccode2) +
                strata(numlib1) +
                frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
                data = data,
                method = c("efron"),
                model = TRUE)

summary(lib1)

lib110 <- coxph(formula = Surv(lib110_start, lib110_end, liberalization_110) ~
                  kin +
                  std_polyarchy_lag1+
                  majorannintra_contig_lag1+
                  ihs_aidgdp_5ma_wins+
                  ihs_gdppc_lag1 +
                  neggdpshock_lag1+
                  ihs_pop_lag1 +
                  cwbroad_lag1 +
                  ihs_iterate_lag1 +
                  ihs_tradeshare_lag1 +
                  
                  cluster(ccode2) +
                  strata(numlib110) +
                  frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                  frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
               data = data,
               method = c("efron"),
               model = TRUE)

summary(lib110)

lib120 <- coxph(formula = Surv(lib12_start, lib12_end, liberalization_12) ~
                  kin +
                  std_polyarchy_lag1+
                  majorannintra_contig_lag1+
                  ihs_aidgdp_5ma_wins+
                  ihs_gdppc_lag1 +
                  neggdpshock_lag1+
                  ihs_pop_lag1 +
                  cwbroad_lag1 +
                  ihs_iterate_lag1 +
                  ihs_tradeshare_lag1 +
                  
                  cluster(ccode2) +
                  strata(numlib12) +
                  frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                  frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
               data = data,
               method = c("efron"),
               model = TRUE)

summary(lib120)

lib130 <- coxph(formula = Surv(lib13_start, lib13_end, liberalization_13) ~
                  kin +
                  std_polyarchy_lag1+
                  majorannintra_contig_lag1+
                  ihs_aidgdp_5ma_wins+
                  ihs_gdppc_lag1 +
                  neggdpshock_lag1+
                  ihs_pop_lag1 +
                  cwbroad_lag1 +
                  ihs_iterate_lag1 +
                  ihs_tradeshare_lag1 +
                  
                  cluster(ccode2) +
                  strata(numlib13) +
                  frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                  frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
                data = data,
                method = c("efron"),
                model = TRUE)

summary(lib130)

lib140 <- coxph(formula = Surv(lib14_start, lib14_end, liberalization_14) ~
                  kin +
                  std_polyarchy_lag1+
                  majorannintra_contig_lag1+
                  ihs_aidgdp_5ma_wins+
                  ihs_gdppc_lag1 +
                  neggdpshock_lag1+
                  ihs_pop_lag1 +
                  cwbroad_lag1 +
                  ihs_iterate_lag1 +
                  ihs_tradeshare_lag1 +
                  
                  cluster(ccode2) +
                  strata(numlib14) +
                  frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                  frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
                data = data,
                method = c("efron"),
                model = TRUE)

summary(lib140)

lib150 <- coxph(formula = Surv(lib15_start, lib15_end, liberalization_15) ~
                  kin +
                  std_polyarchy_lag1+
                  majorannintra_contig_lag1+
                  ihs_aidgdp_5ma_wins+
                  ihs_gdppc_lag1 +
                  neggdpshock_lag1+
                  ihs_pop_lag1 +
                  cwbroad_lag1 +
                  ihs_iterate_lag1 +
                  ihs_tradeshare_lag1 +
                  
                  cluster(ccode2) +
                  strata(numlib15) +
                  frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                  frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
                data = data,
                method = c("efron"),
                model = TRUE)

summary(lib150)

libplot <- multiplot(lib50, lib60, lib70, lib80, lib90, ewlib1, lib1, lib110, lib120, lib130, lib140, lib150, secret.weapon=T, innerCI=2, coefficients = c("kin"), sort = c("magnitude"), names=c("+0.50 SD Liberalization", "+0.60 SD Liberalization", "+0.70 SD Liberalization", "+0.80 SD Liberalization", "+0.90 SD Liberalization", "+1 SD Liberalization \n EW Index", "+1 SD Liberalization \n ICW Index", "+1.10 SD Liberalization", "+1.20 SD Liberalization", "+1.30 SD Liberalization", "+1.40 SD Liberalization", "+1.50 SD Liberalization"), horizontal = F, color = c("black", "black" , "black" , "black" , "black"), zeroColor ="red", textAngle = 0, title= "", xlab = "Standardized Coefficient  \n", ylab="")
libplot

data$kin <- NULL

## Figure A.28

m1b <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               exc1inc2_reg_lag1 +
               std_polyarchy_lag1+
               majorannintra_contig_lag1+
               ihs_aidgdp_5ma_wins+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               policy_top+
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
             frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m1b)

m2b <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               exc1inc2_reg_lag1 +
               std_polyarchy_lag1+
               majorannintra_contig_lag1+
               ihs_aidgdp_5ma_wins+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               british +french +portugeuse +belgian +italian + russian+
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m2b)

m3b <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               exc1inc2_reg_lag1 +
               std_polyarchy_lag1+
               majorannintra_contig_lag1+
               ihs_aidgdp_5ma_wins+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               ihs_unemployilo_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m3b)

m4b <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               exc1inc2_reg_lag1 +
               std_polyarchy_lag1+
               majorannintra_contig_lag1+
               ihs_aidgdp_5ma_wins+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               ihs_natresgdp_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m4b)

m5b <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               exc1inc2_reg_lag1 +
               std_polyarchy_lag1+
               majorannintra_contig_lag1+
               ihs_aidgdp_5ma_wins+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               chinaproj_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m5b)

m6b <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               exc1inc2_reg_lag1 +
               std_polyarchy_lag1+
               majorannintra_contig_lag1+
               ihs_aidgdp_5ma_wins+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               rightgov_lag1 + leftgov_lag1 + nationalistgov_lag1+
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m6b)

m7b <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
                exc1inc2_reg_lag1 +
               std_polyarchy_lag1+
               majorannintra_contig_lag1+
               ihs_aidgdp_5ma_wins+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
                stratdisp_1+
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m7b)

m8b <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               exc1inc2_reg_lag1 +
               std_polyarchy_lag1+
               majorannintra_contig_lag1+
               ihs_aidgdp_5ma_wins+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               coldwar+
               post911+
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m8b)

m9b <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               exc1inc2_reg_lag1 +
               std_polyarchy_lag1+
               majorannintra_contig_lag1+
               ihs_aidgdp_5ma_wins+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               numdomlaws+
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m9b)

m10b <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               exc1inc2_reg_lag1 +
                std_polyarchy_lag1+
                majorannintra_contig_lag1+
                ihs_aidgdp_5ma_wins+
                ihs_gdppc_lag1 +
                neggdpshock_lag1+
                ihs_pop_lag1 +
                cwbroad_lag1 +
                ihs_iterate_lag1 +
                ihs_tradeshare_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
             data = data.law,
             method = c("efron"),
             model = TRUE)

summary(m10b)

m11b <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               exc1inc2_reg_lag1 +
                std_polyarchy_lag1+
                majorannintra_contig_lag1+
                ihs_aidgdp_5ma_wins+
                ihs_gdppc_lag1 +
                neggdpshock_lag1+
                ihs_pop_lag1 +
                cwbroad_lag1 +
                ihs_iterate_lag1 +
                ihs_tradeshare_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
             data = data.1990,
             method = c("efron"),
             model = TRUE)

summary(m11b)

m12b <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               exc1inc2_reg_lag1 +
                std_polyarchy_lag1+
                majorannintra_contig_lag1+
                ihs_aidgdp_5ma_wins+
                ihs_gdppc_lag1 +
                neggdpshock_lag1+
                ihs_pop_lag1 +
                cwbroad_lag1 +
                ihs_iterate_lag1 +
                ihs_tradeshare_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
             data = data.indp10,
             method = c("efron"),
             model = TRUE)

summary(m12b)

controlplot <- multiplot(m1b, m2b, m3b, m4b, m5b, m6b, m7b, m8b, m9b, m10b, m11b, m12b, secret.weapon=T,  innerCI=2, coefficients = c("exc1inc2_reg_lag1"), sort = c("magnitude"), names=c("Control 1: Policy Level", "Control 2: Colonial History", "Control 3: Unemployment", "Control 4: Natural Resource Rents", "Control 5: Chinese Aid", "Control 6: Ruling Party Orientation", "Control 7: Strategic Displacement", "Control 8: Cold War and Post-9/11", "Control 9: # of Policies", "Subset 1: At Least 1 Policy", "Subset 2: Post-1990 Subset", "Subset 3: Post-Independence \n (10+ Years) Subset"), horizontal = F, color = c("black", "black" , "black" , "black" , "black"), zeroColor ="red", textAngle = 0, title= "", xlab = "Standardized Coefficient  \n", ylab="")
controlplot

data$kin <- NULL

## FIGURE A.29

data$kin <- data$exc1inc2_reg_lag1

km1500 <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               kin +
                 std_polyarchy_lag1+
                 majorannintra_contig_lag1+
                 ihs_aidgdp_5ma_wins+
                 ihs_gdppc_lag1 +
                 neggdpshock_lag1+
                 ihs_pop_lag1 +
                 cwbroad_lag1 +
                 ihs_iterate_lag1 +
                 ihs_tradeshare_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(km1500)

data$kin <- NULL
data$kin <- data$exc1inc2_950_lag1

km950 <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
                 kin +
                 std_polyarchy_lag1+
                 majorannintra_contig_lag1+
                 ihs_aidgdp_5ma_wins+
                 ihs_gdppc_lag1 +
                 neggdpshock_lag1+
                 ihs_pop_lag1 +
                 cwbroad_lag1 +
                 ihs_iterate_lag1 +
                 ihs_tradeshare_lag1 +
                 
                 cluster(ccode2) +
                 strata(numlib1) +
                 frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                 frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
               data = data,
               method = c("efron"),
               model = TRUE)

summary(km950)

data$kin <- NULL
data$kin <- data$exc1inc2_reg2000_lag1

km2000 <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
                  kin +
                  std_polyarchy_lag1+
                  majorannintra_contig_lag1+
                  ihs_aidgdp_5ma_wins+
                  ihs_gdppc_lag1 +
                  neggdpshock_lag1+
                  ihs_pop_lag1 +
                  cwbroad_lag1 +
                  ihs_iterate_lag1 +
                  ihs_tradeshare_lag1 +
                  
                  cluster(ccode2) +
                  strata(numlib1) +
                  frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                  frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
                data = data,
                method = c("efron"),
                model = TRUE)

summary(km2000)

data$kin <- NULL
data$kin <- data$exc1inc2_3000_lag1

km3000 <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
                  kin +
                  std_polyarchy_lag1+
                  majorannintra_contig_lag1+
                  ihs_aidgdp_5ma_wins+
                  ihs_gdppc_lag1 +
                  neggdpshock_lag1+
                  ihs_pop_lag1 +
                  cwbroad_lag1 +
                  ihs_iterate_lag1 +
                  ihs_tradeshare_lag1 +
                  
                  cluster(ccode2) +
                  strata(numlib1) +
                  frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                  frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
                data = data,
                method = c("efron"),
                model = TRUE)

summary(km3000)

data$kin <- NULL
data$kin <- data$exc1inc2_reg4000_lag1

km4000 <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
                  kin +
                  std_polyarchy_lag1+
                  majorannintra_contig_lag1+
                  ihs_aidgdp_5ma_wins+
                  ihs_gdppc_lag1 +
                  neggdpshock_lag1+
                  ihs_pop_lag1 +
                  cwbroad_lag1 +
                  ihs_iterate_lag1 +
                  ihs_tradeshare_lag1 +
                  
                  cluster(ccode2) +
                  strata(numlib1) +
                  frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                  frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
                data = data,
                method = c("efron"),
                model = TRUE)

summary(km4000)

data$kin <- NULL
data$kin <- data$exc1inc2_5000_lag1

km5000 <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
                  kin +
                  std_polyarchy_lag1+
                  majorannintra_contig_lag1+
                  ihs_aidgdp_5ma_wins+
                  ihs_gdppc_lag1 +
                  neggdpshock_lag1+
                  ihs_pop_lag1 +
                  cwbroad_lag1 +
                  ihs_iterate_lag1 +
                  ihs_tradeshare_lag1 +
                  
                  cluster(ccode2) +
                  strata(numlib1) +
                  frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                  frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
                data = data,
                method = c("efron"),
                model = TRUE)

summary(km5000)

data$kin <- NULL
data$kin <- data$exc1inc2_reg6000_lag1

km6000 <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
                  kin +
                  std_polyarchy_lag1+
                  majorannintra_contig_lag1+
                  ihs_aidgdp_5ma_wins+
                  ihs_gdppc_lag1 +
                  neggdpshock_lag1+
                  ihs_pop_lag1 +
                  cwbroad_lag1 +
                  ihs_iterate_lag1 +
                  ihs_tradeshare_lag1 +
                  
                  cluster(ccode2) +
                  strata(numlib1) +
                  frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                  frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
                data = data,
                method = c("efron"),
                model = TRUE)

summary(km6000)

data$kin <- NULL
data$kin <- data$exc1inc2_reg7000_lag1

km7000 <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
                  kin +
                  std_polyarchy_lag1+
                  majorannintra_contig_lag1+
                  ihs_aidgdp_5ma_wins+
                  ihs_gdppc_lag1 +
                  neggdpshock_lag1+
                  ihs_pop_lag1 +
                  cwbroad_lag1 +
                  ihs_iterate_lag1 +
                  ihs_tradeshare_lag1 +
                  
                  cluster(ccode2) +
                  strata(numlib1) +
                  frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                  frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
                data = data,
                method = c("efron"),
                model = TRUE)

summary(km7000)

data$kin <- NULL
data$kin <- data$exc1inc2_reg8000_lag1

km8000 <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
                  kin +
                  std_polyarchy_lag1+
                  majorannintra_contig_lag1+
                  ihs_aidgdp_5ma_wins+
                  ihs_gdppc_lag1 +
                  neggdpshock_lag1+
                  ihs_pop_lag1 +
                  cwbroad_lag1 +
                  ihs_iterate_lag1 +
                  ihs_tradeshare_lag1 +
                  
                  cluster(ccode2) +
                  strata(numlib1) +
                  frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                  frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
                data = data,
                method = c("efron"),
                model = TRUE)

summary(km8000)

distanceplot <- multiplot(km950, km1500, km2000, km3000, km4000, km5000, km6000, km7000, km8000, secret.weapon=T, innerCI=2, coefficients = c("kin"), sort = c("magnitude"), names=c("0950 km.", "1500 km.", "2000 km." , "3000 km.", "4000 km.", "5000 km.", "6000 km." , "7000 km.", "8000 km."), horizontal = F, color = c("black", "black" , "black" , "black" , "black"), zeroColor ="red", textAngle = 0, title= "", xlab = "Standardized Coefficient  \n", ylab="")
distanceplot
data$kin <- NULL

## FIGURE A.30

algeria <- subset(data, country!="Algeria")
armenia <- subset(data, country!="Armenia")
azerbaijan <- subset(data, country!="Azerbaijan")
burundi <- subset(data, country!="Burundi")
cameroon <- subset(data, country!="Cameroon")
car <- subset(data, country!="Central African Republic")
djibouti <- subset(data, country!="Djibouti")
iran <- subset(data, country!="Iran")
kazakhstan <- subset(data, country!="Kazakhstan")
kyrgyzstan <- subset(data, country!="Kyrgyzstan")
madagascar <- subset(data, country!="Madagascar")
morocco <- subset(data, country!="Morocco")
senegal <- subset(data, country!="Senegal")
sierraleone <- subset(data, country!="Sierra Leone")
somalia <- subset(data, country!="Somalia")
ssudan <- subset(data, country!="South Sudan")
turkey <- subset(data, country!="Turkey")
uganda <- subset(data, country!="Uganda")
zambia <- subset(data, country!="Zambia")

algeria <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
                   exc1inc2_reg_lag1 +
                   std_polyarchy_lag1+
                   majorannintra_contig_lag1+
                   ihs_aidgdp_5ma_wins+
                   ihs_gdppc_lag1 +
                   neggdpshock_lag1+
                   ihs_pop_lag1 +
                   cwbroad_lag1 +
                   ihs_iterate_lag1 +
                   ihs_tradeshare_lag1 +
                   
                   cluster(ccode2) +
                   strata(numlib1) +
                   frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                   frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
                 data = algeria,
                 method = c("efron"),
                 model = TRUE)

summary(algeria)

armenia <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
                   exc1inc2_reg_lag1 +
                   std_polyarchy_lag1+
                   majorannintra_contig_lag1+
                   ihs_aidgdp_5ma_wins+
                   ihs_gdppc_lag1 +
                   neggdpshock_lag1+
                   ihs_pop_lag1 +
                   cwbroad_lag1 +
                   ihs_iterate_lag1 +
                   ihs_tradeshare_lag1 +
                   
                   cluster(ccode2) +
                   strata(numlib1) +
                   frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                   frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
                 data = armenia,
                 method = c("efron"),
                 model = TRUE)

summary(armenia)

azerbaijan <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
                      exc1inc2_reg_lag1 +
                      std_polyarchy_lag1+
                      majorannintra_contig_lag1+
                      ihs_aidgdp_5ma_wins+
                      ihs_gdppc_lag1 +
                      neggdpshock_lag1+
                      ihs_pop_lag1 +
                      cwbroad_lag1 +
                      ihs_iterate_lag1 +
                      ihs_tradeshare_lag1 +
                      
                      cluster(ccode2) +
                      strata(numlib1) +
                      frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                      frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
                    data = azerbaijan,
                    method = c("efron"),
                    model = TRUE)

summary(azerbaijan)

burundi <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
                   exc1inc2_reg_lag1 +
                   std_polyarchy_lag1+
                   majorannintra_contig_lag1+
                   ihs_aidgdp_5ma_wins+
                   ihs_gdppc_lag1 +
                   neggdpshock_lag1+
                   ihs_pop_lag1 +
                   cwbroad_lag1 +
                   ihs_iterate_lag1 +
                   ihs_tradeshare_lag1 +
                   
                   cluster(ccode2) +
                   strata(numlib1) +
                   frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                   frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
                 data = burundi,
                 method = c("efron"),
                 model = TRUE)

summary(burundi)

cameroon <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
                    exc1inc2_reg_lag1 +
                    std_polyarchy_lag1+
                    majorannintra_contig_lag1+
                    ihs_aidgdp_5ma_wins+
                    ihs_gdppc_lag1 +
                    neggdpshock_lag1+
                    ihs_pop_lag1 +
                    cwbroad_lag1 +
                    ihs_iterate_lag1 +
                    ihs_tradeshare_lag1 +
                    
                    cluster(ccode2) +
                    strata(numlib1) +
                    frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                    frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
                  data = cameroon,
                  method = c("efron"),
                  model = TRUE)

summary(cameroon)

car <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               exc1inc2_reg_lag1 +
               std_polyarchy_lag1+
               majorannintra_contig_lag1+
               ihs_aidgdp_5ma_wins+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
             data = car,
             method = c("efron"),
             model = TRUE)

summary(car)

djibouti <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
                    exc1inc2_reg_lag1 +
                    std_polyarchy_lag1+
                    majorannintra_contig_lag1+
                    ihs_aidgdp_5ma_wins+
                    ihs_gdppc_lag1 +
                    neggdpshock_lag1+
                    ihs_pop_lag1 +
                    cwbroad_lag1 +
                    ihs_iterate_lag1 +
                    ihs_tradeshare_lag1 +
                    
                    cluster(ccode2) +
                    strata(numlib1) +
                    frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                    frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
                  data = djibouti,
                  method = c("efron"),
                  model = TRUE)

summary(djibouti)

iran <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
                exc1inc2_reg_lag1 +
                std_polyarchy_lag1+
                majorannintra_contig_lag1+
                ihs_aidgdp_5ma_wins+
                ihs_gdppc_lag1 +
                neggdpshock_lag1+
                ihs_pop_lag1 +
                cwbroad_lag1 +
                ihs_iterate_lag1 +
                ihs_tradeshare_lag1 +
                
                cluster(ccode2) +
                strata(numlib1) +
                frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
              data = iran,
              method = c("efron"),
              model = TRUE)

summary(iran)

kazakhstan <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
                      exc1inc2_reg_lag1 +
                      std_polyarchy_lag1+
                      majorannintra_contig_lag1+
                      ihs_aidgdp_5ma_wins+
                      ihs_gdppc_lag1 +
                      neggdpshock_lag1+
                      ihs_pop_lag1 +
                      cwbroad_lag1 +
                      ihs_iterate_lag1 +
                      ihs_tradeshare_lag1 +
                      
                      cluster(ccode2) +
                      strata(numlib1) +
                      frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                      frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
                    data = kazakhstan,
                    method = c("efron"),
                    model = TRUE)

summary(kazakhstan)

kyrgyzstan <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
                      exc1inc2_reg_lag1 +
                      std_polyarchy_lag1+
                      majorannintra_contig_lag1+
                      ihs_aidgdp_5ma_wins+
                      ihs_gdppc_lag1 +
                      neggdpshock_lag1+
                      ihs_pop_lag1 +
                      cwbroad_lag1 +
                      ihs_iterate_lag1 +
                      ihs_tradeshare_lag1 +
                      
                      cluster(ccode2) +
                      strata(numlib1) +
                      frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                      frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
                    data = kyrgyzstan,
                    method = c("efron"),
                    model = TRUE)

summary(kyrgyzstan)

madagascar <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
                      exc1inc2_reg_lag1 +
                      std_polyarchy_lag1+
                      majorannintra_contig_lag1+
                      ihs_aidgdp_5ma_wins+
                      ihs_gdppc_lag1 +
                      neggdpshock_lag1+
                      ihs_pop_lag1 +
                      cwbroad_lag1 +
                      ihs_iterate_lag1 +
                      ihs_tradeshare_lag1 +
                      
                      cluster(ccode2) +
                      strata(numlib1) +
                      frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                      frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
                    data = madagascar,
                    method = c("efron"),
                    model = TRUE)

summary(madagascar)

morocco <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
                   exc1inc2_reg_lag1 +
                   std_polyarchy_lag1+
                   majorannintra_contig_lag1+
                   ihs_aidgdp_5ma_wins+
                   ihs_gdppc_lag1 +
                   neggdpshock_lag1+
                   ihs_pop_lag1 +
                   cwbroad_lag1 +
                   ihs_iterate_lag1 +
                   ihs_tradeshare_lag1 +
                   
                   cluster(ccode2) +
                   strata(numlib1) +
                   frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                   frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
                 data = morocco,
                 method = c("efron"),
                 model = TRUE)

summary(morocco)

senegal <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
                   exc1inc2_reg_lag1 +
                   std_polyarchy_lag1+
                   majorannintra_contig_lag1+
                   ihs_aidgdp_5ma_wins+
                   ihs_gdppc_lag1 +
                   neggdpshock_lag1+
                   ihs_pop_lag1 +
                   cwbroad_lag1 +
                   ihs_iterate_lag1 +
                   ihs_tradeshare_lag1 +
                   
                   cluster(ccode2) +
                   strata(numlib1) +
                   frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                   frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
                 data = senegal,
                 method = c("efron"),
                 model = TRUE)

summary(senegal)

sierraleone <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
                       exc1inc2_reg_lag1 +
                       std_polyarchy_lag1+
                       majorannintra_contig_lag1+
                       ihs_aidgdp_5ma_wins+
                       ihs_gdppc_lag1 +
                       neggdpshock_lag1+
                       ihs_pop_lag1 +
                       cwbroad_lag1 +
                       ihs_iterate_lag1 +
                       ihs_tradeshare_lag1 +
                       
                       cluster(ccode2) +
                       strata(numlib1) +
                       frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                       frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
                     data = sierraleone,
                     method = c("efron"),
                     model = TRUE)

summary(sierraleone)

somalia <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
                   exc1inc2_reg_lag1 +
                   std_polyarchy_lag1+
                   majorannintra_contig_lag1+
                   ihs_aidgdp_5ma_wins+
                   ihs_gdppc_lag1 +
                   neggdpshock_lag1+
                   ihs_pop_lag1 +
                   cwbroad_lag1 +
                   ihs_iterate_lag1 +
                   ihs_tradeshare_lag1 +
                   
                   cluster(ccode2) +
                   strata(numlib1) +
                   frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                   frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
                 data = somalia,
                 method = c("efron"),
                 model = TRUE)

summary(somalia)

ssudan <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
                  exc1inc2_reg_lag1 +
                  std_polyarchy_lag1+
                  majorannintra_contig_lag1+
                  ihs_aidgdp_5ma_wins+
                  ihs_gdppc_lag1 +
                  neggdpshock_lag1+
                  ihs_pop_lag1 +
                  cwbroad_lag1 +
                  ihs_iterate_lag1 +
                  ihs_tradeshare_lag1 +
                  
                  cluster(ccode2) +
                  strata(numlib1) +
                  frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                  frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
                data = ssudan,
                method = c("efron"),
                model = TRUE)

summary(ssudan)

turkey <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
                  exc1inc2_reg_lag1 +
                  std_polyarchy_lag1+
                  majorannintra_contig_lag1+
                  ihs_aidgdp_5ma_wins+
                  ihs_gdppc_lag1 +
                  neggdpshock_lag1+
                  ihs_pop_lag1 +
                  cwbroad_lag1 +
                  ihs_iterate_lag1 +
                  ihs_tradeshare_lag1 +
                  
                  cluster(ccode2) +
                  strata(numlib1) +
                  frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                  frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
                data = turkey,
                method = c("efron"),
                model = TRUE)

summary(turkey)

uganda <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
                  exc1inc2_reg_lag1 +
                  std_polyarchy_lag1+
                  majorannintra_contig_lag1+
                  ihs_aidgdp_5ma_wins+
                  ihs_gdppc_lag1 +
                  neggdpshock_lag1+
                  ihs_pop_lag1 +
                  cwbroad_lag1 +
                  ihs_iterate_lag1 +
                  ihs_tradeshare_lag1 +
                  
                  cluster(ccode2) +
                  strata(numlib1) +
                  frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                  frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
                data = uganda,
                method = c("efron"),
                model = TRUE)

summary(uganda)

zambia <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
                  exc1inc2_reg_lag1 +
                  std_polyarchy_lag1+
                  majorannintra_contig_lag1+
                  ihs_aidgdp_5ma_wins+
                  ihs_gdppc_lag1 +
                  neggdpshock_lag1+
                  ihs_pop_lag1 +
                  cwbroad_lag1 +
                  ihs_iterate_lag1 +
                  ihs_tradeshare_lag1 +
                  
                  cluster(ccode2) +
                  strata(numlib1) +
                  frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                  frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
                data = zambia,
                method = c("efron"),
                model = TRUE)

summary(zambia)

outlierplot <- multiplot(algeria, armenia, azerbaijan, burundi, cameroon, car, djibouti, iran, kazakhstan, kyrgyzstan, madagascar, morocco, senegal, sierraleone, somalia, ssudan, turkey, uganda, zambia, secret.weapon=T, innerCI=2, coefficients = c("exc1inc2_reg_lag1"), sort = c("magnitude"), names=c("Algeria", "Armenia", "Azerbaijan" , "Burundi", "Cameroon", "Central African Republic", "Djibouti" , "Iran", "Kazakhstan", "Kyrgyzstan", "Madagascar", "Morocco", "Senegal", "Sierra Leone", "Somalia", "South Sudan", "Turkey", "Uganda", "Zambia"), horizontal = F, color = c("black", "black" , "black" , "black" , "black"), zeroColor ="red", textAngle = 0, title= "", xlab = "Standardized Coefficient  \n", ylab="")
outlierplot
